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^ ■ We consider the problem of the existence of a dynamical barrier of "mass" that needs to be 

O ' excited on a lattice site to lead to the formation and subsequent persistence of localized modes for a 

, nonlinear Schrodinger lattice. We contrast the existence of a dynamical barrier with its absence in 

the static theory of localized modes in one spatial dimension. We suggest an energetic criterion that 
provides a sufficient, but not necessary, condition on the amplitude of a single-site initial condition 
required to form a solitary wave. We show that this effect is not one-dimensional by considering 
its two-dimensional analog. The existence of a sufficient condition for the excitation of localized 
modes in the non-integrable, discrete, nonlinear Schrodinger equation is compared to the dynamics 
of excitations in the integrable, both discrete and continuum, version of the nonlinear Schrodinger 
equation. 
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I. INTRODUCTION 



In the past few years, there has been an explosion of interest in discrete models that has been summarized in 
K*" I numerous recent reviews This growth has been spurted by numerous apphcations of dynamical lattice nonlinear 
' models in areas as diverse as the nonlinear optics of waveguide arrays 0] , the dynamics of Bose-Einstein condensates 
in periodic potentials [2], micro-mechanical models of cantilever arrays or even simple models of the complex 
dynamics of the DNA double strand [1]. Perhaps the most prototypical model among those that emerge in these 
settings is the, so-called, discrete nonlinear Schrodinger equation (DNLS) @. DNLS may arise as a direct model, as 
a tight binding approximation, or even as an envelope wave expansion: it is, arguably, one of the most ubiquitous 
models in the nonlinear physics of dispersive, discrete systems. 
I In at least one of these settings (namely, in the nonlinear optics of waveguide arrays with the focusing nonlinearity), 
the feature that will be of interest to the present work has been observed experimentally. In particular, it has been 
noted, to the best of our knowledge firstly in Ref. Q, that when an injected beam of light into one waveguide had 
low intensity, then the beam dispersed through quasi-linear propagation. On the other hand, in the same work, 
$^ ' experiments with high intensity of the input beam led to the first example of formation of discrete solitary waves in 
5^ , waveguide arrays. A very similar "crossover" from linear to nonlinear behavior was also observed very recently in arrays 
of waveguides with the defocusing nonlinearity Q . The common feature of both works is that they used the DNLS 
equation as the supporting model to illustrate this behavior at a theoretical/numerical level. However, this crossover 
phenomenon is certainly not purely discrete in nature. Perhaps the most famous example of a nonlinear wave equation 
that possesses such a threshold is the integrable continuum nonlinear Schrodinger equation 9] . Specifically, it is well- 
known that, e.g., in the case of a square barrier of initial conditions of amplitude Vq and width L, the product VqL 
determines the nature of the resulting soliton, and if it is sufficiently small the initial condition disperses without the 
formation of a solitonic structure [To| . On the other hand, the existence of the threshold is not a purely one-dimensional 
feature either. For instance, experiments on the formation of solitary waves in two-dimensional photorefractive crystals 
show that low intensities lead to diffraction, whereas higher intensities induce localization pUfl^. Moreover, similar 
phenomena were observed even in the formation of higher-order excited structures such as vortices (as can be inferred 
by carefully inspecting the results of Refs. [H, HI])- It should be mentioned that the latter field of light propagation 
in photorefractive crystals is another major direction of current research in nonlinear optics; see, e.g., Ref. [15] for a 
recent review. 

This crossover behavior between linear and nonlinear dynamics may be understood qualitatively rather simply. In 
the case of power law nonlinearities of order p, which are relevant in these settings, a small intensity 0(e), where 



2 



e <C 1, yields a nonlinear contribution 0{eP) that is negligible with respect to the linear terms of the equation. On 
the other hand, if e S> 1, the opposite will be true and the nonlinear terms will dominate the linear ones, yielding 
essentially nonlinear behavior. A key question regards the details of this crossover and what determines its more 
precise location for an appropriately parametrized initial condition. This is the question we address herein. We argue 
that the problem related to the experiments described above reduces, at the mathematical level, to a DNLS equation 
with a Kronecker-(5 initial condition parametrized by its amplitude. Then, a well defined value of the initial-state 
amplitude exists such that initial states with higher amplitude always give rise to localized modes. The condition 
may be determined by comparing the energy of the initial state with the energy of the localized excitations that the 
model supports. This sufficient, but not necessary, condition for the formation of localized solitary waves provides an 
intuitively and physically appealing interpretation of the dynamics that is in very good agreement with our numerical 
observations. We also consider variants of this process in different settings: for reasons of completeness, we present it 
also in the continuum NLS equation, noting the significant differences that the latter case has from the present one. 
As yet another example of very different (from both its non-integrable sibling and its continuum limit) dynamical 
behavior, we also present the case of the integrable discrete NLS (so-called Ablowitz-Ladik model. In 

addition to the one-dimensional DNLS lattice, we also consider the two-dimensional case where the role of both 
energy and beam power (mathematically the squared Z^-norm) become apparent. We should note here that our tool 
of choice for visualizing the "relaxational process" (albeit in a Hamiltonian system) of the initial condition will be 
energy-power diagrams. Such diagrams have proven very helpful in visualizing the dynamics of initial conditions in a 
diverse host of nonlinear wave equations. In particular, they have been used in the nonlinear homogeneous systems 
such as birefringent media and nonlinear couplers as is discussed in Chapters 7 and 8 of Ref. [l^ . They have also 
been used in a form closely related to the present work (bu t in the continuum case; see also the discussion below) 
for general nonlinearities in dispersive wave equations in [20|, while they have been used to examine the migration of 
localized excitations in DNLS equations in [2l| . 

Our presentation is structured as follows. In section II, we present the analytically tractable theory of the integrable 
"relatives" of the present model: we review the known theory for the NLS model and develop its analog for the 
integrable discrete NLS case. Then, in section III, we present our analytical and numerical results in the one- and 
two-dimensional DNLS equation. In the last section, we summarize our findings and present our conclusions, as well 
as highlight some important questions for future studies. 

II. THRESHOLD CONDITIONS FOR THE INTEGRABLE NLS MODELS 
A. The Continuum NLS Model 

For reasons of completeness of the presentation and to compare and contrast the results of the non-integrable case 
that is at the focus of the present work, we start by summarizing the threshold conditions for the continuum NLS 
model [lOj . For the focusing NLS equation 

iut = -^u^x - \ufu, (1) 

with squared barrier initial data 

-(-'0)-{o,' tth-er^Se^ (2) 

(the inverse of) the transmission coefficient, Sii{E), which is the first entry of the scattering matrix, is given by 

Sii{E) = v{E) cos{2v{E)L) - iE sin{2v{E)L), (3) 

with v{E) — y^EP + V^ where E is the spectral parameter and Vq the amplitude of the barrier. It is well-known 
that the number of zeros of this coefficient represents the number of solitons produced by the square barrier initial 
condition p^ . 

It can be proved that the roots of this equation are purely imaginary. (This initial condition satisfies the single-lobe 
conditions of Klaus-Shaw potentials, from which it follows that the eigenvalues are purely imaginary ^25|). Let us 
define 77 > and use E = irj. Then, Eq. ([3]) becomes 

- cos + V sin = 0. (4) 
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We can verify that Eq. (|4|) does not have any roots (i.e., leads to no sohtons in Eq. Ilj) if 



V, < \. (5) 



Furthermore, the condition to generate n sohtons, i.e., so that Eq. (jl]) has n roots is 



(2n-l)| <2yoi< (27i + l)|, (6) 

or, equivalently, the count of eigenvalues is given by 

-VoL-\<n<-VaL+\. (7) 

The limit Vq — > cxd together with L — )■ can be reached if we impose 2V{)L = const. In this instance, the number of 
eigenvalues stays the same. 

B. The Ablowitz-Ladik Model 

We now turn to the integrable discretization of Eq. ([T]) and examine its dynamics. The one-dimensional integrable, 
discrete nonlinear Schrodinger model (so-called Ablowitz-Ladik (AL) model [IB, [13, [3 ) reads: 

iUn = -^(w«+l + -"ri-l - 2u„) - ^|u„P (u„+i -|- U„_i) . (8) 

In this case, there exists a Lax pair of linear operators [1^, [13, [ill 

£„ - Z + M„, (9) 
Bn = (^-^) D + ]^{ZMn- Z-^Mn-i) -]^DMnM,,^i, (10) 
with the definitions for the matrices 



where z is the spectral parameter, and Un = Un{t) is a solution of the equation. These two operators ^ and (flO)) 
define the system of differential-difference equations 

^'„+i = (12) 

i^-^n = Bn-^n. (13) 

dr 

for a complex matrix function Then, the compatibility condition of Eqs. (I12p and (I13p 



m—n+l 



(i.e., the Lax equation) becomes the AL model. [/„ = Un{t) is referred to as the potential of the AL eigenvalue 
problem. 

For Un decaying rapidly at ±oo, and for n — > ±oo, from Eq. (|12p we have: 



We normalize this type of solutions as follows: let denote the solution of Eq. ([72]) such that 

~ Z" as n—> +0G, 

and let $„ be the solution of Eq. ([12]) such that 

$„ Z" as n ^ —CO. 
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'^n and are known as the Jost functions. Each of these forms a system of hnearly independent solutions of the 
AL eigenvalue problem p^ . These sets of solutions are inter-related by the scattering matrix S{z), 

$„ = *„5(z). (14) 

The first column of this equation is given by 

($i)„ = 5n(z)(*i)„ + 52i(^)(*2)„, (15) 

where (<&i)n denotes the first column of Similar definitions apply to ('l'i)n and (^'2)n- Since $„ ~ Z" as n ^ — oo, 
then 

(cl>i)„^z"(^j 

To obtain decay, ($i)n when n — > — oo, we demand 

1^1 > 1- (16) 

On the other hand, 

/z" \ 

= ((«'l)«,(*2)rO ~ ( n-^oo. 

Therefore, if \z\ > 1 then (5'i)n ^ oo and (^'2)n ^ 0, as 7i oo. Now, from Eq. it follows that (<i>i)n — !■ oo as 
n ^ oo, unless S'ii(z) = 0. 

We therefore seek solutions zi, Z2, . . . , z^r of the equation 

^ii(zfe) = 0, fc = l,2,...,7V, (17) 

such that |zfc| > 1. Then, 

{^l)nizk) ^ S2lizk){^2)n{zk), fc = 1, 2, . . . , TV. 

From this, it follows that ($i)„(zfe) decays at ±00: 

(<i>i)„(zi:) — > as n ^ ±00. 

It is then said that ($i)n(zfc) is an eigenfunction (fc = 1, 2, . . . , N), with corresponding eigenvalue z^. 
In the case of C/„(t = 0) = f/o(5„,„o, the Jost function is 

= Z"-i(Z + Mo)$o, for n>l, (18) 

with Mo 

Furthermore, = for ?i > 1 and $0 is the identity matrix: $0 = I- Hence, Eq. reads: 

= ^'„Z"i(Z + A/o), for n > 1. 
Its comparison with Eq. ([H]) leads to a scattering matrix 

S{z) = Z-i(Z + Afo). 
We thus obtain that the transmission coefficient 

^11 (^) = 1, 

which never vanishes. This means that the one-site potential (i.e., a single-site initial condition) does not admit 
solitonic solutions, independently of the amplitude Uq of initial excitation. This theoretical result has also been 
confirmed by numerical simulations for different values of Uq, always leading to dispersion of the solution (not shown 
here). 
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III. THRESHOLD CONDITIONS FOR THE NON-INTEGRABLE DNLS MODEL 



We now turn to the non-integrable DNLS lattice that has the general form: 

-eA2W„ - 



|2cr„ 



(20) 



where is a complex field (corresponding to the envelope of the electric field in optics, or the mean field wavefunction 
in optical lattice wells in the BECs), A2U„ = (un+i + Un-i — 2u„) is the discrete Laplacian, and e is the ratio of the 
tunneling strength to the nonlinearity strength. Using a scaling invariance of the equation, we can scale e = 1, through 
rescaling t ^ et and m„ Unj Importantly for our considerations, this model has an underlying Hamiltonian 
(which is also a conserved quantity in the dynamics) of the form: 
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(21) 



The only other known conserved quantity in the dynamics of Eq. (j20p is the beam power (in Bose- Einstein condensates, 
the normalized number of atoms) P = 'Ylm l^^nP- 

The specific problem that we examine here addresses the following question. Consider a "compactum" of mass 
Un ~ ASn,o', what is the critical value of A that is necessary for this single-site initial condition to excite a localized 
mode? That such a threshold definitely exists is illustrated in Fig. [TJ The case of the leftmost panel is subcritical, 
leading to the discrete dispersion of the initial datum. This follows the well-known t~^^'^ amplitude decay which is 
implied by the solution of the problem in the absence of the nonlinearity: 



Unit) = Az"J„(2t), 



(22) 



where Jn is a Bessel function of order n. Notice that this can also be shown to be consistent with the findings of 
[2^ 1. On the other hand, the rightmost panel shows a nonlinearity-dominated regime with the rapid formation of a 
solitary wave strongly localized around n = 0. In the intermediate case of the middle panel, the system exhibits a 
long oscillatory behavior reminiscent of a separatrix between the basins of attraction of the two different regimes. 
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FIG. 1: Subcritical case (left), critical (middle), and supercritical (right) initial single-site excitations on the lattice. In all 
cases, the initial condition is Un = A5n,o, with A = 1 in the left, 2 in the middle, and 2.5 in the right panels. In each case the 
top panel shows the space-time contour plot of the evolution. The bottom panel shows the dynamical evolution of |it„|^ for 
sites n = 0, 1, 2 (solid, dashed and dash-dotted lines, respectively). For the leftmost case the inset shows the same evolution in 
a log-log plot and a decay for comparison. Clearly, the damped oscillation of the field modulus has an envelope of t~^^^. 



We argue herein that the essence of this separatrix lies in the examination of the stationary (localized) states of 
the model. Such standing wave solutions of the form m„ — exp(iAi)w„, which are exponentially decaying for Vn as a 
function of n, can be found for arbitrary frequency A (and arbitrary power P in one spatial dimension). This is a 
well-known result in one dimension; see e.g., Ref. [2j|. On the other hand, the single-site initial condition discussed 
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above has an energy of: 

Hss = 2A^ - (23) 
cr + i 

where the subscript denotes smgle site. Figure [2] summarizes succinctly the power dependence of the energy for these 
two cases for cr = 1. Both the energy of the stationary solutions as a function of their power and the single-site 
energy as a function of single-site power {Pss = A'^), are shown. Note that the two curves Hss{Pss) and H{P) do not 
intersect (except at the trivial point H = P = 0) since for e ^ single-site states are not stationary ones. 

The examination of Fig. [2] provides information on the existence of a sufficient condition for the formation of a 
localized mode and on the dynamics of the single-site initial condition. Figure [2] shows that localized solutions exist 
for arbitrarily small values of the input power P, and that the energy of the localized states is negative. This implies 
that the crucial quantity to determine the fate of the process is the energy H and not the power P. The role of the 
power will become more evident in the two-dimensional setting mentioned below. Moreover, if the system starts at 
a given point on the curve defined by Hgs = 2Pss — Pg^^ /{a + 1), due to conservation of total H and total P, it 
can only end up in a stationary state in the quadrant H < H^s and P < Pgs- That is to say, some of the initial 
energy and power are typically "shed off" in the form of radiation (i.e., converted to other degrees of freedom which 
is the only way that "effective dissipation" can arise in a purely Hamiltonian system), so that the initial condition 
can "relax" to the pertinent final configuration. As mentioned above, a localized solution with the same power as 
that of the initial condition exists for arbitrary A. However, emergence of a localized mode occurs only for those 
initial conditions whose core energy (i.e., the energy of a region around the initially excited state) is negative, after 
the profile is "reshaped" by radiating away both energy and power. 

Therefore, if Hss < 0, then the compactum of initial data will always yield a localized excitation: this inequality 
provides the sufficient condition for the excitation of solitary waves. The condition on the energy, in turn, provides 
a condition on the single-site amplitude that leads to the formation of solitary waves, namely, solitary waves always 
form if A> A* with 

A* = [2 {a + 1)]^ . (24) 

For the case of ct = 1 considered in Figs. [Hand [21 this amplitude value is A* — 2 in agreement with our numerical 
observations of Fig. [TJ Whether an initial state with Hss > yields a localized state depends on the explicit system 
dynamics, corroborating the observation that the previous energetic condition is a sufficient, but not necessary, 
condition. 
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FIG. 2: The blue solid line shows the energy H versus the power P of the discrete solitary wave solutions. Above it, the green 
solid line shows the energy versus power of the initial condition, obtained from Eq. (|23|l and Pss ~ A^. The horizontal line 
denotes H — and its intersection with the initial condition curve defines the single-site, initial amplitude A* . The dynamical 
evolution of three different supercritical initial states with A — 2.1, A = 2.3, and A — 2.5 is shown by dashed (black) lines; see 
also discussion in the text. 

One should make a few important observations here. Firstly, we note the stark contrast of the non-integrable 
discrete model and both of its integrable (continuum and discrete) counterparts. In the (singular) continuum limit, it 
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is possible to excite a single soliton or a multi-soliton depending on the barrier height and width. On the other hand, 
in the integrable discretization one-site excitation never leads to solitary wave formation, contrary to what is the case 
here where either none or one solitary wave may arise, depending on the amplitude of the initial one-site excitation. 

Restricting our consideration to the non-integrable model, we observe that even though a localized solution with 
the same power as that of the initial condition exists for arbitrary A, formation of a localized mode will always occur 
only for A > A* , as determined by the previous energetic criterion. Secondly, the answer that the formation always 
occurs for A > A* generates two interconnected questions: what is the threshold for the formation of the localized 
mode, and given an initial A> A* , which one among the mono-parametric family of solutions will the dynamics of the 
model select as the end state of the system ? [It should be noted in connection to the latter question that single-site 
initial conditions have always been found in our numerical simulations to give rise to at most a single-site-centered 
solution i.e., multipulses cannot be produced by this process.] Some examples of this dynamical process are illustrated 
in Fig. [51 where the energy and power of a few sites (typically 20-40) around the originally excited one are measured 
as a function of time and are parametrically plotted in the H — P plane. As it should, the relevant curve starts 
from the Hss — Pss curve, and asymptotically approaches, as a result of the dynamical evolution, the H — P curve 
of the stationary states of the system. However, the relaxation process happens neither at fixed energy, nor at fixed 
power. Instead, it proceeds through a more complex, dynamically selected pathway of loss of both H and P to relax 
eventually to one of the relevant stationary states. This is shown for three different values of super-critical amplitude 
in Fig. = 2.1, 2.3 and 2.5). We have noticed (numerically) that the loss of energy and power, at least in the initial 
stages of the evolution happens at roughly the same rate, resulting in dH/dP « const. However, the later stages of 
relaxation no longer preserve this constant slope. It is worthwhile to highlight here that the fact that the process of 
formation of a localized mode is neither equienergetic, nor does it occur at fixed power is something that has been 
previously observed in the continuum version of the system in (20| (cf. with the discussion in p. 6095 therein). 

The dynamics of the system in the H — P space can be qualitatively understood by considering the frequencies 
A associated with the instantaneous H, P values along the system trajectories. At each instantaneous energy H 
a frequency A^ may be defined as the frequency of a single-site breather stationary state with energy H. Such 
a frequency is unique as the stationary-state energy is a monotonically decreasing function of the frequency [2^ . 
Similarly, a frequency Ap may be uniquely associated with the instantaneous power P, i.e., Ap is the frequency of a 
stationary state with power P. Note, however, that the stationary-state power is a monotonically increasing function 
of A. At the final stationary state where the system relaxes Ah = Ap = —dH/dP. Hence, the system trajectories 
are such that Ah increases (consequently the energy decreases) and Ap decreases (the power decreases). The final 
stationary state is reached when the two frequencies become equal, the point where they meet depending on their rate 
of change along the trajectory, i.e., on their corresponding "speeds" along the trajectory. Nevertheless, the precise 
mechanism of selection of the particular end state (i.e., of the particular "equilibrium A") that a given initial state 
will result in remains a formidable outstanding question that would be especially interesting to address in the future. 

We now turn to the two-dimensional variant of the above one-dimensional non-integrable lattice. Equation (j20p 
remains the same, but for the two-dimensional field Un,m, and the discrete Laplacian becomes the five-point stencil 
A2W„,m = {un+i,m + -I- -l- w„,m-i " ^Un,m)- For the initial condition u„^m = A5„^o^m,o, the Hgs - Pss 

curve is given by: 

i/,, -4F,,-P;/V(a + 1). (25) 

In d-dimensions the first term would be 2dPss. We once again find the branch of standing waves of the equation. 
While this branch of solutions is also well-known psj , there are some important differences with the one-dimensional 
case. Firstly, a stable and an unstable branch of solutions exists (per the well-known Vakhitov-Kolokolov criterion); in 
the wedge-like curve indicating the standing waves only the lower energy branch is stable. Furthermore, the maximal 
energy of the solutions is no longer H = 0, but finite and positive (in fact, for tr = 1, it is Hq « 1.85). Finally, 
solutions no longer exist for arbitrarily low powers, but they may only exist above a certain power (often referred to 
as the excitation threshold [l^ll^, HH). We can now appreciate the impact of these additional features in the right 
panel of Fig. [S] Comparing Hss with Hq we find two solutions: one with A 0.7 and one with A « 2.73. However, 
for the lower one there are no standing wave excitations with the corresponding power (this reveals the role of the 
power in the higher dimensional problem). Hence, the relevant amplitude that determines the sufficient condition 
in the two-dimensional case is the latter. Indeed, the case with A = 2.65 shown in the left panel of Fig. [3] shows 
spatio-temporal diffraction, while that of A = 2.75 in the middle panel illustrates robust localization, indicating the 
super-critical nature of the latter case. Notice that the pathway of the dynamics is presented in the right panel for 
three super-critical values of yl = 2.75, A — 2.9 and A = 3.1. Once again, the dynamics commences on the Hss ~ Pss 
curve, as it should, eventually relaxing on the stable standing wave curve. As before, the initial dynamics follows a 
roughly constant dH/dP, but the relaxation becomes more complex at later phases of the evolution. 
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FIG. 3: The left panel shows the spatio-temporal evolution of \u„,m{t)\'^ for the sub-critical case with amplitude A = 2.65. 
The middle panel shows the weakly super-critical case of A = 2.75. Finally the right panel is analogous to Fig. (2] showing the 
H ~ P diagram for the solution branch (wedge- like blue solid line), the Hss — Pss graph of the initial conditions (green solid 
line) and the trajectories of three super-critical cases in H ~ P space for A = 2.75, 2.9 and 3.1 (dashed black lines). The red 
horizontal line represents the maximal energy for which solutions are found to exist, namely Ho ~ 1.85. 

IV. CONCLUSIONS AND FUTURE CHALLENGES 

The above study has examined the presence of a sharp crossover between the linear and nonlinear dynamics of a 
prototypical dynamical lattice model such as the discrete nonlinear Schrodinger equation. This crossover has already 
been observed in media with the focusing nonlinearity ^7,] (as considered here). It has also been observed very recently 
in media with the defocusing nonlinearity Q. The latter can be transformed into the former under the so-called 
staggering transformation u„ = (— where Wn is the field in the defocusing case. As a result, the solitary waves 
of the defocusing problem discussed in Q will be "staggered" (i.e., of alternating phase between neighboring sites), yet 
the phenomenology discussed above will persist. The crossover was quantified on the basis of an energetic comparison 
of the initial-state energy with the branch of corresponding stable localized solutions "available" in the model. A 
sufficient, but not necessary, condition for the excitation of a localized mode based on the initial-state, single-site 
amplitude was discussed, and it was successfully tested in numerical simulations. Similar findings were obtained in 
the two-dimensional analog of the problem: a crossover behaviour dictated by the energy was found, but the crossover 
was also affected by the power and its excitation thresholds. Furthermore, these results were contrasted with the case 
of the continuum version of the model, where depending on the strength of the excitation, also multi-solitons can be 
obtained and with the integrable discrete Ablowitz-Ladik model where no single-site excitation can produce a solitary 
wave, independently of the excitation amplitude. 

However, a number of interesting questions emerge from these findings that are pertinent to future studies. Perhaps 
the foremost among them concerns how the dynamics "selects" among the available steady-state excitations with 
energy and power below that of the initial condition the one to which the dynamical evolution leads. According to 
our results, this evolution is neither equienergetic, nor power-preserving (see also [2(1]), hence it would be extremely 
interesting to identify the leading physical principle which dictates it. Another question generalizing the more principal 
one asked herein concerns the excitation of multiple sites (possibly at the same amplitude), starting with two, and 
inquiring the nature of the resulting state of the system. On a related note, perhaps this last question is most suitable 
to be addressed first in the context of the integrable model where the formulation presented herein can be generalized 
and yield definitive answers for the potential excitation of solitons. Studies along these directions are currently in 
progress and will be reported in future publications. 
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